The goal of this analysis is to explore patterns in the distribution and relative abundance of krill from Beth Phillips’ work. Then, we will spatially match it to our eulachon eDNA samples in order to use krill abundacne as a predictor in eulachon SDMs.
Import the krill backscatter data. The data represent relative krill densities between 50-300 m in the water column, along all transects that collected acoustic data during the years 2019, 2021, and 2023. In the data, ID is a combination of Transect number and Interval (horizontal 0.5 nmi bin, with numbering starting at 1 at the start of each transect), and the min/max depth of each 10m bin. All of the zeroes are still included, along with positive NASC values.
kdf <- read_csv(here('data','HakeSurvey_KrillNASC_0.5nmix10m_cells_2019-2023.csv'))
## Rows: 470392 Columns: 12
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): ID, Species, Transect
## dbl (8): Year, Interval, Layer_depth_min, Layer_depth_max, NASC, Lat, Lon, ...
## date (1): Date
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
Import the eulachon samples
d_obs <- read_rds(here('Data','eulachon qPCR 2019 and 2021 samples clean.rds')) %>%
mutate(Ct=replace_na(Ct,0)) %>% # delta-models in sdmTMB need this
mutate(utm.lon.km=utm.lon.m/1000,
utm.lat.km=utm.lat.m/1000)
Background map
#background map
pred.crs <- terra::rast(here('Data','raster_grid_blake','fivekm_grid.tif')) %>% st_crs()
coast <- ne_states(country='United States of America',returnclass = 'sf') %>%
filter(name %in% c('California','Oregon','Washington','Nevada')) %>%
st_transform(crs = pred.crs)
Bathymetry
# bathy
limits.for.map <- d_obs %>%
st_as_sf(coords=c('lon','lat'),crs=4326) %>%
st_bbox()+c(-1,-1,1,1) #add an extra degree on each side for buffer
b = getNOAA.bathy(lon1 = limits.for.map["xmin"],
lon2 = limits.for.map[ "xmax" ],
lat1 = limits.for.map["ymin"],
lat2 = limits.for.map["ymax"],
resolution = 1,keep = TRUE)
## File already exists ; loading 'marmap_coord_-127.429775;33.4438333333333;-119.62;49.4493333333333_res_1.csv'
# make into a dataframe of contours for ggplotting
bdf <- fortify(b) %>%
contoureR::getContourLines(levels=c(0, -100, -200, -500, -1000,-1500)) %>%
st_as_sf(coords=c("x","y"),crs=4326) %>%
st_transform(pred.crs) %>%
mutate(x=st_coordinates(.)[,1],y=st_coordinates(.)[,2]) %>%
st_set_geometry(NULL)
ksf <- kdf %>% st_as_sf(coords=c("Lon","Lat"),crs=4326) %>% st_transform(pred.crs)
# raw
ksf %>%
ggplot(aes(NASC))+
geom_density()+theme_minimal()
# natural log
ksf %>%
ggplot(aes(log(NASC)))+
geom_density()+theme_minimal()+
labs(x="log (NASC)")
## Warning: Removed 434303 rows containing non-finite values (`stat_density()`).
# by depth bin
ksf %>%
mutate(big_depth_category=cut(Layer_depth_min,10)) %>%
ggplot(aes(log(NASC),color=factor(Year)))+
geom_density(size=1.5)+theme_minimal()+
scale_color_aaas()+
labs(x="log (NASC)",title="Distribution by Depth Category")+
facet_wrap(~big_depth_category)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## Removed 434303 rows containing non-finite values (`stat_density()`).
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
log (NASC) looks pretty normal, with a heavy left tail (zeroes). When
faceted by depth bin, it’s clear that in general, there were more krill
in 2019 than in 2021 or 2023.
latdepth_plot <- kdf %>%
mutate(Layer_depth_mid=Layer_depth_min+5) %>%
mutate(NASC=na_if(NASC,0)) %>%
ggplot(aes(Layer_depth_mid,Lat,size=log(NASC+1),color=log(NASC+1)))+
geom_point()+
scale_color_viridis(option="magma",na.value='gray80')+
scale_size_area(max_size=4,na.value=0.5)+
facet_wrap(~Year)+
theme_minimal()+theme(panel.background = element_rect(color="black",fill=NA))
latdepth_plot
This combined figure of latitude, depth, and NASC shows the change in
depth distribution by latitude across years. In it, it is clear that
there are differences in not only the depth distribution of krill among
years (as shown in the last two plots), but also there are hotspots in
different places along the coast (e.g., the presence of krill at 40N in
2019 but not really in either subsequent year; the large surface
hotspots that occur in 2021 but not in other years at ~38N and ~47N,
etc.)
Similar to the above, we can parse the NASC outputs by transect to look at patterns.
krill_by_x <- kdf %>%
group_by(Transect) %>%
mutate(xLatMean=mean(Lat)) %>%
arrange(xLatMean) %>%
group_by(Transect,xLatMean) %>%
nest() %>%
mutate(p=purrr::map2(Transect,data,function(x,df){
df %>%
ggplot(aes(Layer_depth_min,NASC,color=factor(Year),shape=factor(Year)))+
geom_point()+
labs(color="Year",shape="Year",x="Depth Layer",y="NASC",title=x)+
theme_minimal()+
# geom_smooth(se=F)+
theme(panel.border = element_rect(fill=NA,color='black'))+
scale_color_manual(breaks=c("2019","2021","2023"),values=viridis_pal(option="G",end=0.8)(3),drop=F)+
scale_shape_manual(breaks=c("2019","2021","2023"),values=c(16,17,18),drop=F)
}))
Latitudinal layout of transects (for reference)
krill_by_x %>% ggplot(aes(fct_reorder(Transect,xLatMean),xLatMean))+geom_point()+
theme_minimal()+
theme(axis.text.x=element_text(angle=90,vjust=0))+
ggtitle("Mean Transect Latitude")
Plots of NASC by depth across transects, arranged South to North by
average latitude of each transect
plts <- krill_by_x %>% pull(p)
plot_grid(plotlist=plts[1:20],nrow=4)
plot_grid(plotlist=plts[21:40],nrow=4)
plot_grid(plotlist=plts[41:60],nrow=4)
plot_grid(plotlist=plts[61:80],nrow=4)
plot_grid(plotlist=plts[81:100],nrow=4)
plot_grid(plotlist=plts[101:120],nrow=4)
plot_grid(plotlist=plts[121:134],nrow=4)
Some interesting interannual changes are obvious here, both depth shifts between years and magnitude of the signal.
Similar to the plots above, these are coastwide maps of krill NASC by lumped depth category. In the plots, larger and “hotter” (yellow/red) points indicate higher NASC values. Many of the same patterns (e.g., greater abundance in 2019 in many places) are apparent.
NASC_by_depth_map <- ksf %>%
mutate(xs=str_remove_all(Transect,"x") %>% as.numeric()) %>%
mutate(big_depth_category=cut(Layer_depth_min,10)) %>%
mutate(NASC=na_if(NASC,0)) %>%
filter(xs<95) %>%
arrange(big_depth_category) %>%
group_by(big_depth_category) %>%
nest() %>%
mutate(p=purrr::map2(big_depth_category,data,function(d,df){
df %>%
ggplot()+
geom_sf(aes(color=log(NASC+1),size=log(NASC+1)),shape=1)+
facet_wrap(~Year,nrow=1)+
scale_color_viridis(option="C",limits=c(0,11),na.value='gray80')+
scale_size_area(limits=c(0,11),max_size=8,na.value=0.5)+
theme_minimal()+
ggtitle(d)+
theme(panel.border=element_rect(color='black',fill=NA))
}))
NASC_by_depth_map$p
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
With these comparisons among years in hand, we can spatially match the NASC data to the eulachon eDNA samples. Do spatial matching because the transect numbers are seemingly all mixed up.
kdf <- kdf %>% mutate(transect_num=str_remove_all(Transect,"x") %>% as.integer())
matchkey <- d_obs %>% distinct(date,year,month,day,station,lat,lon,depth_cat) %>%
mutate(transect1=str_split_i(station,"-",1) %>% as.integer(),
transect2=str_split_i(station,"_",1) %>% str_remove_all("x") %>% as.integer()) %>%
mutate(interval1=str_split_i(station,"-",2) %>% as.integer(),
interval2=str_split_i(station,"_",2) %>% as.integer()) %>%
mutate(transect_match=coalesce(transect1,transect2),
interval_match=coalesce(interval1,interval2)) %>%
left_join(kdf,by=c("date"="Date","year"="Year","month"="Month","depth_cat"="Layer_depth_min","transect_match"="transect_num","interval_match"="Interval"))
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `transect1 = str_split_i(station, "-", 1) %>% as.integer()`.
## Caused by warning in `str_split_i(station, "-", 1) %>% as.integer()`:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `interval1 = str_split_i(station, "-", 2) %>% as.integer()`.
## Caused by warning in `str_split_i(station, "-", 2) %>% as.integer()`:
## ! NAs introduced by coercion
Matching by transect is not working well. Let’s make a flat (depth-integrated) version and then match by location
kdf_2d <- kdf %>%
group_by(Year,ID,Lat,Lon) %>%
summarise(sumNASC=sum(NASC)) %>%
ungroup() %>%
mutate(latbin=cut(Lat,5))
kdf_2d_maps <- kdf_2d %>%
mutate(sumNASC=na_if(sumNASC,0)) %>%
unite(facet,Year,latbin) %>%
group_by(facet) %>%
nest() %>%
mutate(p=purrr::map2(facet,data,function(l,d){
d %>%
ggplot(aes(Lon,Lat,color=log(sumNASC),size=log(sumNASC)))+
geom_point(shape=1)+
scale_size(range=c(0,8),limits=c(0,10))+
scale_color_viridis(option="C",limits=c(0,10),na.value="gray80")+
labs(title=facet)+
coord_equal()+
theme_minimal()+
theme(panel.border = element_rect(color='black',fill=NA))
}))
kdf_2d_maps$p
## [[1]]
## Warning: Removed 780 rows containing missing values (`geom_point()`).
##
## [[2]]
## Warning: Removed 1240 rows containing missing values (`geom_point()`).
##
## [[3]]
## Warning: Removed 1013 rows containing missing values (`geom_point()`).
##
## [[4]]
## Warning: Removed 1321 rows containing missing values (`geom_point()`).
##
## [[5]]
## Warning: Removed 1026 rows containing missing values (`geom_point()`).
##
## [[6]]
## Warning: Removed 550 rows containing missing values (`geom_point()`).
##
## [[7]]
## Warning: Removed 1217 rows containing missing values (`geom_point()`).
##
## [[8]]
## Warning: Removed 641 rows containing missing values (`geom_point()`).
##
## [[9]]
## Warning: Removed 1001 rows containing missing values (`geom_point()`).
##
## [[10]]
## Warning: Removed 858 rows containing missing values (`geom_point()`).
##
## [[11]]
## Warning: Removed 1011 rows containing missing values (`geom_point()`).
##
## [[12]]
## Warning: Removed 1173 rows containing missing values (`geom_point()`).
##
## [[13]]
## Warning: Removed 891 rows containing missing values (`geom_point()`).
##
## [[14]]
## Warning: Removed 1101 rows containing missing values (`geom_point()`).
##
## [[15]]
## Warning: Removed 1103 rows containing missing values (`geom_point()`).
These plots are hard to parse visually—still trying to think of the best
way to display them.
Let’s match eDNA and krill by year and location. There are multiple ways we can conceptualize that this matching might work, depending on how much smoothing or averaging we want to do, across both latitude/longitude and depth. For example, we can match each eDNA sample to its nearest neighbor in lat/lon space, then extract the krill NASC value from the closest depth to our eDNA sample. Alternatively, we could do the same, but use a depth-integrated version of krill, as we calculated in the previous step. Finally, we could average krill NASC among, e.g., the five nearest neighbors to each eDNA sample, or among all samples within some distance of the eDNA sample.
In order to stay flexible to the above options, we’ll write a flexible function to find matches for individual eDNA samples.
# krill data locations, converted to a flat projection
krill_locs <- kdf_2d %>%
st_as_sf(coords=c("Lon","Lat"),crs=4326,remove=F) %>%
st_transform(pred.crs) %>%
mutate(utm.lon.m=st_coordinates(.)[,1],utm.lat.m=st_coordinates(.)[,2])%>%
st_set_geometry(NULL)
# this function takes a row/observation from the eDNA dataset, and finds its nn nearest neighbors on the same transect
match_krill <- function(dfrow,nn=5){
# transect number
x_num <- dfrow$transect_num
# location of eDNA sample
this_loc <- dfrow %>% dplyr::select(utm.lon.m,utm.lat.m)
# krill data to query (same transect, same year)
krill_to_compare <- kdf %>% filter(Year==dfrow$year,transect_num==x_num)%>%
st_as_sf(coords=c("Lon","Lat"),crs=4326) %>% st_transform(pred.crs)
# locations of the krill data
krill_locs <- krill_to_compare %>%
distinct(ID,.keep_all = T)
# determine the nn nearest neighbors
nns <- nn2(st_coordinates(krill_locs),this_loc,k=nn)
# find the krill data that were identified as neighbors
krill_matched <- krill_locs %>%
st_set_geometry(NULL) %>%
slice(as.numeric(nns$nn.idx)) %>%
# include the distance (here, in meters)
mutate(nn.dist=nns$nn.dists %>% as.numeric) %>%
# squared/inverse distance weighting (for later use)
mutate(dist.weight=1/(nn.dist^2)) %>%
dplyr::select(Year,ID,nn.dist,dist.weight)
# finally, use the matched locations to filter the krill data, then return
krill_final <- krill_to_compare %>%
st_set_geometry(NULL) %>%
left_join(krill_matched,by=join_by(Year,ID)) %>%
filter(!is.na(nn.dist)) %>%
dplyr::select(Year,ID,Layer_depth_min,Layer_depth_max,NASC,nn.dist,dist.weight)
krill_final
}
# test the above function
# edna_krill_matched_5nn <- d_obs %>%
# slice(6700:6725) %>%
# mutate(transect1=str_split_i(station,"-",1) %>% as.integer(),
# transect2=str_split_i(station,"_",1) %>% str_remove_all("x") %>% as.integer()) %>%
# mutate(transect_num=coalesce(transect1,transect2)) %>%
# dplyr::select(-transect1,-transect2) %>%
# mutate(rn=row_number()) %>%
# group_by(rn) %>%
# nest() %>%
# mutate(krill=purrr::map(data,match_krill)) %>%
# ungroup()
# test <- edna_krill_matched_5nn %>% unnest(cols=c(data))
# Run on entire dataset
library(tictoc)
tic("Matching krill nearest neighbors")
edna_krill_matched_5nn <- d_obs %>%
mutate(transect1=str_split_i(station,"-",1) %>% as.integer(),
transect2=str_split_i(station,"_",1) %>% str_remove_all("x") %>% as.integer()) %>%
mutate(transect_num=coalesce(transect1,transect2)) %>%
dplyr::select(-transect1,-transect2) %>%
mutate(rn=row_number()) %>%
group_by(rn) %>%
nest() %>%
mutate(krill=purrr::map(data,match_krill)) %>%
ungroup()
toc()
# this join takes ~30-35m
We end up with a nested data frame of krill data (5 nearest neighbors on the transect, and all depth slices) assigned to each “observation” (row) of the eDNA data.
write_rds(edna_krill_matched_5nn,here('model output','edna_krill_matched_5nn.rds'))
# if you're not re-running the above, can just import here
edna_krill_matched_5nn <- read_rds(here('model output','edna_krill_matched_5nn.rds'))
Options:
edna_krill_metrics <- edna_krill_matched_5nn %>%
mutate(depth_slice=purrr::map(data,function(df)df %>% pull(depth_cat))) %>%
# depths that will match krill "Layer_depth_min"
mutate(depth_match=case_when(
depth_slice==0 ~ 50L,
depth_slice==25 ~ 50L,
depth_slice==50 ~ 50L,
depth_slice==100 ~ 100L,
depth_slice==150 ~ 150L,
depth_slice==200 ~ 200L,
depth_slice==300 ~ 290L,
depth_slice==500 ~ 290L,
)) %>%
# one more depth match- if the depth slice is deeper than the deepest krill sample for that location, use max krill depth instead
mutate(depth_match=purrr::map2(depth_match,krill,function(depth,df){
ifelse(depth>max(df$Layer_depth_min),max(df$Layer_depth_min),depth)
})) %>%
# first metric- closest depth, nearest neighbor
mutate(k1=purrr::map2_dbl(krill,depth_match,function(df,depth){
out <- df %>% filter(Layer_depth_min==depth) %>% filter(nn.dist==min(nn.dist)) %>% pull(NASC) %>% as.numeric()
out
})) %>%
# second metric-closest depth slice, averaged across neighbors with inverse distance weighting
mutate(k2=purrr::map2_dbl(krill,depth_match,function(df,depth){
nndf <- df %>% filter(Layer_depth_min==depth)
out <- weighted.mean(nndf$NASC,w=nndf$dist.weight)
out
})) %>%
# third metric-closest neighbor, integrate across depth
mutate(k3=purrr::map2_dbl(krill,depth_match,function(df,depth){
nndf <- df %>% filter(nn.dist==min(nn.dist))
out <- sum(nndf$NASC)
out
})) %>%
# fourth metric-integrate across depth, then average across neighbors with inverse distance weighting
mutate(k4=purrr::map2_dbl(krill,depth_match,function(df,depth){
nndf <- df %>% group_by(nn.dist,dist.weight) %>% summarise(sumNASC=sum(NASC)) %>% ungroup()
out <- weighted.mean(nndf$sumNASC,w=nndf$dist.weight)
out
})) %>%
# fifth metric-integrate across depth, then average across neighbors with distance weighting (but not squared distances)
mutate(k5=purrr::map2_dbl(krill,depth_match,function(df,depth){
nndf <- df %>% group_by(nn.dist,dist.weight) %>% summarise(sumNASC=sum(NASC)) %>% ungroup()
out <- weighted.mean(nndf$sumNASC,w=1/nndf$nn.dist)
out
}))
#unnest
edna_krill_metrics <- edna_krill_metrics %>%
dplyr::select(data,k1:k5) %>%
unnest(cols=c(data))
write_rds(edna_krill_metrics,here('model output','edna_krill_metrics.rds'))
# if you need to read (i.e., you don't need or want to re-run the slow step above)
edna_krill_metrics <- read_rds(here('model output','edna_krill_metrics.rds'))
Another option is to “smooth” the krill data using an SDM. This is a little odd, since we want the krill data to feed into an SDM for eulachon, but nevertheless, let’s calculate it and see what it looks like. We will fit a simple model in sdmTMB with only spatial and spatial-temporal latent factors (i.e., no environmental covariates), then project onto a fine-scale grid. We can also project the model at the exact locations of the eDNA data, thereby adding the model-smoothed krill densities as yet another option for a metric.
We make this model based on the depth-integrated (2D) data.
# re-project 2D krill data to a flat projection
kdf_2d_utm <- kdf_2d %>%
st_as_sf(coords=c("Lon","Lat"),crs=4236) %>%
st_transform(pred.crs) %>%
mutate(x=st_coordinates(.)[,1],y=st_coordinates(.)[,2]) %>%
st_set_geometry(NULL) %>%
mutate(xkm=x/1000,ykm=y/1000)
# make the INLA mesh for sdmTMB
mesh <- make_mesh(kdf_2d_utm, c("xkm", "ykm"), cutoff = 20)
plot(mesh)
## Fit
Now we fit a simple model in sdmTMB:
kfit <- sdmTMB(
sumNASC ~ 0+as.factor(Year),
data = kdf_2d_utm,
mesh = mesh,
time="Year",
spatiotemporal = "IID",
family = tweedie(link='log')
)
#AIC 157419 for sumNASC ~ 1
# AIC 157420.5 for sumNASC ~0+as.factor(Year)
The model converged with no warnings. We can look at the model fit and basic parameter estimation
kfit
## Spatiotemporal model fit by ML ['sdmTMB']
## Formula: sumNASC ~ 0 + as.factor(Year)
## Mesh: mesh (isotropic covariance)
## Time column: Year
## Data: kdf_2d_utm
## Family: tweedie(link = 'log')
##
## coef.est coef.se
## as.factor(Year)2019 1.74 0.39
## as.factor(Year)2021 1.67 0.41
## as.factor(Year)2023 1.09 0.40
##
## Dispersion parameter: 15.41
## Tweedie p: 1.72
## Matérn range: 39.91
## Spatial SD: 2.98
## Spatiotemporal IID SD: 3.92
## ML criterion at convergence: 78702.255
##
## See ?tidy.sdmTMB to extract these values as a data frame.
tidy(kfit,effects = 'ran_pars')
## # A tibble: 5 × 3
## term estimate std.error
## <chr> <dbl> <dbl>
## 1 range 39.9 3.23
## 2 phi 15.4 0.142
## 3 sigma_O 2.98 0.214
## 4 sigma_E 3.92 0.149
## 5 tweedie_p 1.72 0.00215
# model residuals
kfit_resids <- residuals(kfit)
qqnorm(kfit_resids)
qqline(kfit_resids)
We can project this krill model onto a spatial grid and map it.
grid.pred <- read_rds(here('Data','prediction_grid_5km_sdmTMB.rds')) %>%
rename(xkm=utm.lon.km,ykm=utm.lat.km) %>%
filter(bathy.bottom.depth>-500) %>%
mutate(Year=2019) %>%
replicate_df(time_name="Year",time_values=c(2019,2021,2023))
kfit.preds <- predict(kfit,newdata=grid.pred)
Here’s a map of predictions
make_map <- function(dat, column) {
ggplot(dat, aes(xkm, ykm, fill = {{ column }})) +
geom_raster() +
coord_fixed()+
theme_minimal()+
theme(panel.border = element_rect(color='black',fill=NA))
}
kpred.map <- make_map(kfit.preds,est)+
scale_fill_viridis(option="C")+
facet_wrap(~Year)+
labs(title="Krill Predictions",fill="log (NASC)")
kpred.map
Map of fixed effects only. This will look silly because the only fixed effect we have is the effect of year.
kpred.rf.map <- make_map(kfit.preds,est_non_rf)+
scale_fill_viridis(option="C")+
facet_wrap(~Year)+
labs(title="Krill Predictions (fixed effects only)",fill="log (NASC)")
kpred.rf.map
Map of spatial random effects only
kpred.rf.map <- make_map(kfit.preds,omega_s)+
scale_fill_viridis(option="C")+
labs(title="Krill Predictions (spatial random effects only)",fill="log (NASC)")
kpred.rf.map
Map of spatiotemporal effects only (i.e., this is how the consistent spatial effects above vary across years)
kpred.st.map <- make_map(kfit.preds,epsilon_st)+
scale_fill_viridis(option="C")+
facet_wrap(~Year)+
labs(title="Krill Predictions (spatiotemporal random effects only)",fill="log (NASC)")
kpred.st.map
# And an index of total krill abundance...
# area=25 because these are 5km grid cells
kfit.preds.obj <- predict(kfit,newdata=grid.pred,return_tmb_object = T)
tic("Making Abundance Index")
kidx <- get_index(kfit.preds.obj,area=25)
toc()
kidx %>%
ggplot(aes(Year, log_est, ymin=log_est,ymax=upr))+
geom_ribbon(fill='gray50')+geom_line()
We can match these modeled outputs to the eulachon data as well, and see how related the output is to our other metrics.
# Predict at locations of eDNA data
kfit.preds.eul.pts <- d_obs %>%
rename(xkm=utm.lon.km,ykm=utm.lat.km) %>%
mutate(Year=year) %>%
predict(kfit,newdata=.)
edna_krill_metrics$k6 <- exp(kfit.preds.eul.pts$est)
# nearest neighbor matching (deprecated, probably don't need this because we're just predicting at the observation points)
# d_obs_sf <- d_obs %>% st_as_sf(coords=c("utm.lon.m","utm.lat.m"),crs=pred.crs)
# kfit.preds.sf <- kfit.preds %>% st_as_sf(coords=c("x","y"),crs=pred.crs) %>%
# # only need to join the estimate of krill
# dplyr::select(Year,est) %>%
# mutate(k6=exp(est))
#
# d_obs_kfit_join <- purrr::map_df(unique(d_obs_sf$year),function(y){
# d_obs_sf %>% filter(year==y) %>% st_join(kfit.preds.sf %>% filter(Year==y),st_nn,k=1)
# })
#
# edna_krill_metrics$k6 <- d_obs_kfit_join$k6
Are krill and eDNA correlated?
edna_krill_cor <- edna_krill_metrics %>%
filter(Ct>0,depth_cat<300) %>%
dplyr::select(Ct,copies_ul,k1:k6) %>%
# only positive estimates, and only eulachon-relevant samples
cor(use = "complete.obs")
corrplot(edna_krill_cor,method='number')
# pairs plot, including zeroes
edna_krill_pairs_zeroes <- edna_krill_metrics %>%
dplyr::select(Ct,k1:k6) %>%
GGally::ggpairs()+
ggtitle("Pairs plots, all eDNA observations")
edna_krill_pairs_zeroes
# pairs plot, excluding zeroes for Ct
edna_krill_pairs <- edna_krill_metrics %>%
filter(Ct>0, depth_cat<300) %>%
dplyr::select(Ct,k1:k6) %>%
GGally::ggpairs()+
ggtitle("Pairs plots, for all Ct>0")
edna_krill_pairs
# pairs plot, excluding zeroes for Ct AND krill
edna_krill_pairs_positive_only <- edna_krill_metrics %>%
filter(Ct>0, k1>0, depth_cat<300) %>%
dplyr::select(Ct,k1:k6) %>%
GGally::ggpairs()+
ggtitle("Pairs plots, for all Ct>0 AND krill>0")
edna_krill_pairs_positive_only